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Abstract. We study analytically and numerically the stability of the pressure-less, viscously spreading accretion ring. We show 
that the ring is unstable to small non-axisymmetric perturbations. To perform the perturbation analysis of the ring we use a 
stretching transformation of the time coordinate. We find that to 1st order, one-armed spiral structures, and to 2nd order addi- 
tionally two-armed spiral features may appear Furthermore, we identify a dispersion relation determining the instability of the 
ring. The theoretical results are confirmed in several simulations, using two different numerical methods. These computations 
prove independently the existence of a secular spiral instability driven by viscosity, which evolves into persisting leading and 
trailing spiral waves. Our results settle the question whether the spiral structures found in earlier simulations of the spreading 
ring are numerical artifacts or genuine instabilities. 
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In the theory of accretion discs, the idealised problem of a 
viscously spreading, pressure-less ring orbiting a central point 
mass on Keplerian orbits is often used to exemplify the main 
features of an evolving thin accretion disc, i.e. the inward 
transpor t of m a ss and outward tra nsport of angular momentum 
( [Pringle 1981; Frank et al. 2002). With the approximation of 
a small kinematic viscosity v that is independent of the sur- 
face density, an analytic solution exists for this problem, stated 
originally by Liist ( 1952| ) and later by Lynden-Bell & Pringle 
(|1974|). 



in this paper that the non-axisymmetric spiral structures result 
from a genuine physical instability of the problem. 



This axisymmetric analytic solution of the viscous dust ring 
is frequently used as test problem for numerical methods devel- 
oped to simulate accretion discs. The tested algorithms cover 
parti cle methods like smoothed particle hydro dynamics (SPH) 
(e.g. iFlebbe e"tal| |1994 Fpeith & Riffert||l999| ) as well as grid- 
based codes like RH2D (e.g. |Kley||l999P ^. 

However, numerical simulations of the evolving ring often 
show the appearance of additional structures in the disc. In the 
majority of cases, these structures consist of non-axisymmetric 
one-armed spirals, but eventually also narrow concentric rings 
appear. Since their first discovery, it is controversially debated 
whether these structure s are numerical artifacts or mathemat- 
ical instabilities. While Maddison et al. ( 1996) show that for 



the axisymmetric case the concentric rings found in SPH sim- 
ulations are numerical artifacts of the method, we demonstrate 
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Recently, Dgilvie (2001 ) analysed radially extended accre- 
tion discs in a work presenting a general model of eccentric 
accretion discs, of which the viscous ring is a simple special 
case. However, because his main objective was much wider, he 
did not elaborate on the stability properties of the ring. 

The remainder of the paper is organised as follows. In the 
next section, we first summarise the basic hydrodynamic equa- 
tions used for our analysis of the evolving disc, and present for 
reference the analytic solution of the viscously spreading dust 
ring. 

In Sect. H we use a special perturbation technique, the time- 
stretching approach, which allows to obtain solutions of the 
spreading ring at different orders of a small expansion param- 
eter. We derive a linear evolution equation for an eccentricity 
function E. It is shown that to 2nd order, the general solution 
allows the development of one-armed as well as two-armed spi- 
ral features. 

In Sect. Q we perform a local stability analysis of the final 
equation for E and deduce a dispersion relation which indicates 
an instability for some viscosity laws. These results are com- 
pared with numerical simulations we present in Sect. |5] using 
two different numerical methods, one grid-based and the other 
particle based. We find in general very good agreement of both 
numerical methods with the predictions of the local analysis. 
Thus, the simulations confirm the existence of a physical insta- 
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bility in the viscously evolving ring problem. Finally, in || we 
give our conclusions. 

2. The equations 

In accretion physics one assumes that the evolution of the disc 
can be modelled by the hydrodynamic equations. In the present 
case, we envisage a thin accretion disc, i.e. the vertical thick- 
ness of the disc is very small in comparison to the radial ex- 
tend. Therefore, the evolution of the surface density Z of the 
disc can be modelled by the vertically averaged hydrodynamic 
equations which reduces the problem to two dimensions and 
formally corresponds to the 2D-version of the hydrodynamic 
equations. 

2.1. Basic hydrodynamic equations 

Because of the axisymmetry of the ring solution, it is conve- 
nient to work in polar coordinates {R, (p). In the presence of a 
central gravitational field originating by the mass Mc, the set 
of the basic hydrodynamic equations consists of the continuity 
equation 

R 



51 1 diRl-VR) 
5f R dR 



dip 



= 



(1) 



the radial component of the Navier Stokes equation 

,,2 \ 



1 diRvs'^a-Rji) 
^ R dR 



1 divJLaR^) 
R 



' dR 

c 

R 



' R^ 



(2) 



and the azimuthal component of the Navier Stokes equation 



i dt 



+ Vr^ + -ir 



dR 



1 d(Rvj:aR^) 



R dtp R 
1 divJLa-^^ 
R 



1 dp, 
R dip 



) 



-O" Rip 



(3) 



dip R 

Here v^; and denote the radial and azimuthal velocity, ps is 
the vertically averaged pressure, and G is the gravitational con- 
stant. Assuming that the bulk viscosity is vanishing, the viscous 
shear tensor cr has the form 

AdvR 2vr 2 1 dv^ 
3~dR ~ IT ~ 3R~d^ ' 



O-RR - -■ 



dv^ 



C Rp 



1 dvR 



R Rdip 



(4) 
(5) 
(6) 



dR 

2 dvR 4 vr 4 1 dv^ 
'3~dR ^ 3"r 3R~d^ 
and V5 is the vertically averaged kinematic viscosity. 

As we only consider cold (pressure-less) discs, we neglect 
the pressure p^ throughout the rest of this paper. 

2.2. Analytic solution of the viscous ring 

For the time evolution of the viscously spreading ring, we as- 
sume Vj - const. Then, with initial condition 



2^init(-R) - 



M 

IttRo 



6(R-Ro) 



(7) 



at f = with total ring mass M, the surface density of the ring 
evolves according to 

V, . . . M 1 (2x\ I l+j^^ 
(e.g. |:.ynden-Bell & Pringlj |1974[ ), where x = R/Rq, t* 



(8) 



J, 

12vsf/7?y, and L is the modified Bessel function to the order 
of 1/4. 

The azimuthal velocity of the ring is equal to Keplerian ve- 
locity. 



^Kepler - 



R 



= OR 



where we have defined the angular velocity 
niR) = 



R' 

The radial velocity of the ring obeys the relation 

_ 3 d 

^spreading 



(9) 



(10) 



(11) 



Sring dR 

It is important to realize, that solution (||) fulfils the hy- 
drodynamic equations given in 2.1 only by approximation, as- 
suming the kinematic viscosity v, being small compared to the 
specific angular momentum ilR^. 

3. Perturbation analysis 

To study the stability properties of the ring solution (^, we 
perform a perturbation analysis of the hydrodynamic equations 
(§, (|) and (|). 

3.1. Stretching transformation 

According to the idealised problem, we assume that the aver- 
aged kinematic viscosity Vj is very small and depends only on 
R. To take this into account, we replace 

Vj = VsiR) - ev{R) where e = const. 1 . (12) 

Then, two different timescales can be distinguished, the vis- 
cous timescale fvisc ~ R^/vs which determines the radial 
spreading of the ring, and the dynamic timescale fdyn ~ 
R/vp - -sJr^/GMc describing the azimuthal motion of the gas 
in the disc. Because of (p^, the dynamic timescale is much 
shorter than the viscous timescale, fdyn/fvisc ~ Vs/ ^IGM^^R - 
evIiQM^) <K 1. Therefore, in addition to the dynamic time t 
we define a viscous time 



(13) 



In principle, the ring may evolve on the dynamic timescale 
as well as on the viscous timescale. Thus we assume for the 
time dependencies of all hydrodynamic quantities 

m = tit, t) , VRit) = hit, t) , v^it) = v^it, t) . (14) 

This leads to a stretching transformation for the time deriva- 
tives. 



dfit) ^ dfit,T) ^ dfit,T) 

dt dt ^ dr 

where / symbolises any of the hydrodynamic quantities. 



(15) 
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3.2. Approach 

To generalise the approach, we do not restrict the stability anal- 
ysis to a perturbation of the ring solution ^ but we start with 
less determined unperturbed functions. The initial condition of 
the problem shall consist of an arbitrary axisymmetric surface 
density distribution Si„it(/?) rotating around a central mass M^. 
Because self gravity of the disc can be neglected, the initial az- 
imuthal velocity shall be equal to the Keplerian velocity (||), 
V™' - VKepier- Then, the following approach is made for the 
expansion 

1.(1, R,(p) 

VR(t,R,ip) 



I.o(t,T,R) + eI.i(t,T,R,<p) 
+ £^t2(t,T,R,ip)+0(e^) 

VRoit, T,R) + evRi (f , T, R, if) 
+ e^vR2{t,T,R,ip)+0{e^) 

v^o(t, ^) + e v^i (f , r, R, ip) 
+ e^v^2{t,T,R,ip) + 0{e^) 



(16) 
(17) 
(18) 



Here it is assumed that the 0th order quantities are axisymmet- 
ric and do not depend on the azimuthal angle <f, and that the 0th 
order azimuthal velocity evolves only very slowly and is inde- 
pendent of the dynamical time t (both assumptions are justified 
by the initial conditions). 

The approaches (|l6|), (|l7|) and (|l8|) are inserted in ^ (3) 
and (Bb considering the stretching transformation and (0, 
and the terms of the resulting equations are ordered according 
to powers of e. 

There are three main differences of the p resent approac h 
compared to previous analyses like those by Ogilvie ( [2001 ). 
The kinema tic visco s ity is assumed to be a function of radius 
only, while Ogilvie ( 2001 ) considered Vj to be a function of 
surface density. Bot h dynam i c and viscous timescale are taken 
into account, while [Ogilvie ( 2001 ) assumed at the outset that 
non of the variables varies on the dynamic timescale. And the 
expansion is carried out up to second order, so that non-linear 
terms appear in the perturbation. Usually, a linear perturbation 
is performed. 

3.3. Zeroth order results 

In 0(1), the continuity equation reads 

dto 1 diRtoVRo) 



+ — - 



= , 



dt R dR 
the radial component of the Navier Stokes equation reads 



RO 



+ Vro 



(19) 



(20) 



dt ' dR R R^ ' 

and the azimuthal component accordingly 

VRodiRv^^ 
R dR 

From (|T]) follows immediately that the 0th order radial veloc- 
ity has to vanish everywhere. 



vro = Q , 



because the alternative solution of ([2l[), vyo °^ 1 /R, does not 
fulfil the initial condition for the azimuthal velocity. 
Inserting ( p^ into ( [1^ and ( pO[ ) results in 



52o 
dt 

and 



= 



V(pO - ^Keple 



(23) 



(24) 



i.e., the 0th order of the azimuthal velocity is equal to the 
Keplerian velocity, and the 0th order of the surface density does 
not depend on the dynamical time, Zq - ^oiT, R)- 

3.4. First order results 

In 0(e), the continuity equation has the form 



+ ^^r- + -T— + 7^ + TT- 

dt dip dr R dR R dip 



0, 



(25) 



and for the Navier Stokes equation results as radial component 
, (26) 



dvRi „/5vri 
and as azimuthal component 



dv^[ Q 
"df ^ 2 



dip 



Vri + 



d [vto yfR] 



So dR 



= 0. 



(27) 



Taking the derivative of ([27[) with respect to t considering 
( [23| ) and solving for dvR\/dt, and taking the derivative of ( p7[ ) 
with respect to ip and solving for dvRi/dip, and inserting the 
results into (pq) yields 



d^Va,l 



5f2 



H-2Q- 



dtdip 



dip- 



2 + n\i = 



(28) 



The general solution of this equation is a linear superposition 
of 



^R[v.ie'^'«^-^')] 



(29) 



with vyi = v'^i(r, R), and where m and w have to obey the rela- 
tion 

(J - 2mQw + d?-n? - = , 
that is 

(jj - Q.{m + 1) . 

Solving (|7|) for vr\ then results in 
v«i=VR: + 5l[vR,e'('«^-^^)] 
with 

3 d 



(30) 



(31) 



(32) 



vr\ 



and 



[vto V^] 



(22) Vri = 2/^^ -mjv^i 



(33) 



(34) 



4 



R. Speith and W. Kley: Stability of the viscously spreading ring 



Inserting (^) and ( |3^ into gives 



dt dip 



dr R dR 



(35) 



. 

+ im—v^i 
K 



J(m(p - cot) 



The first term in brackets on the right hand side of ( p^ ) and 
the terms in parenthesis in the second line of the right hand 
side of ( p5| ) are secular terms for £1. They have to vanish iden- 
tically, otherwise £1 would increase linearly or, for the latter, 
even quadratically with f or (f, respectively. 
Then, ( p5| ) can be solved for £1, 

ii =1R[S,e'('"^-'^f)] 



with 
2i = 



1 



' dR 



(36) 



(37) 



The conditions for the disappearance of the secular terms 
in parenthesis in ( |3^ are 

dm 
dR 







and 



(38) 



^ =0 
dR 



(39) 



Because of relation ( |3C| ) and Q - f2(/?), this can only be 
achieved by 

co = , (40) 

which leads to - 1 . 

The condition for the disappearance of the secular term in 
brackets in ( |3^ ) is 

dto ^ 1 d(RtoVRi) 
dr 



R 
3d_ 
RdR 



dR 



(41) 



This is the diffusion equation of the surface density for the an- 
alytic solution of the viscous dust ring. With initial condition 

and constant viscosity v, the surface density (||) solves ( |4l| ) 
exactly (considering t* - 12vt/R^). 

That is, the surface density of the analytic solution is equal 
to the 0th order term of the perturbation. The azimuthal velocity 
of the analytic solution is also equal to the according 0th order 
term of the perturbation, i.e., equal to the Keplerian velocity 
(24). And the radial velocity of the analytic solution is equal 
to the axisymmetric part of the according 1st order term of the 
perturbation, i.e., relations dill ) and ( ^3| ) are identical. 

Because of (p^), all first order quantities are independent of 
the dynamic time f . Therefore, the first order velocities take the 
form 



v«i = v«i + nivui e'""^] and 



with 



m — ±1 



(42) 
(43) 
(44) 



Thus, if a non-axisymmetric perturbation appears, to first or- 
der it has to be an one-armed spiral structure. Additionally, it 
follows 



\vri\ 



■2\v. 



(45) 



and vri and v^i obey a phase shift of |. 



Defining a function E - E(t,R) similar to Ogilvie (2001) 

(46) 



with 

vri{t,R) = imRQ.(R)E(T,R) 
yields 



v^i(T,R)^-^RniR)E(T,R) 



(47) 



Then, with (pq), the first order surface density takes the form 
£i(T,/?,^) = !R[Ei(T,/?)e™^] , (48) 
where 



RQm 



d(RI.oVRi) I.QV 



0^'R1 



dR 



-R 



d(i:oE) 
dR 



(49) 



Furthermore, an initial phase i^o can be found such that 
E{t,R) = E(T,R)e''Po with E = l^l. 

3.5. Second order results 

Consider the 0th order results ( ^2| ) and ( |24| ) and the 1 st order 
results (H), (|3j) and Then, in 0{£^) one yields for the 
continuity equation 



d^2 ^d^2 ld{RI.oVR2) ^,)d%2 

+ + 7. + TT— 

dt dip R dR R dip 



(50) 



m d 
~ 2RdR 

\ dR 

^ RdR 



, dCLoE) 



dR 



K Vr\ 



dR 



sin(2mi^ + lipo) 

E d(RtoVRi) 
R dR 

COS(t71(p + ipo) 



for the radial component of the Navier Stokes equation 



So 



dvR2 ^(dvR2 

+ ii I — 2v^2 



dt 



dip 



(51) 



.dE 



■ m£()7?Q — sin(m^ + po) 
dr 



+ -toR^a^E 
2 



3E dE\ 

cos(2m(p + 2po) 

4R dR] ^ ^ 



- niQ. 



1 E dE 
4R ^ M 
E I 1 diR^v%) 



R\R^- dR 

1 dtodivto) 



+ 3R 



dHv^o) 



to dR dR 



dR^ 



1^/ ^(vi^^n 

3dR\ dR 2 



4 . d^E 



sin(m^ + (^o) 



4 1 d 
3RdR 
2 
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Inserting ( p6| ) into ( plj ) results in 



3 dR\R^) 



dvRi _ dvRi 

+ Vr\ 



dT 



dR 



where for some transformations the relation ( P3| ) is used, and 
for the azimuthal component of the Navier Stokes equation 



+ Q 



dv^2 1 , 



1 ^ dE 

- -1.qRQ.— cos(mtp + (fo) 

2 OT 

1 T T dE 
+ -m'LoR^QrE— sin(2m^ + 2(^o) 



Q 



1 . ^2:,, 



5y 



- - -vEo + -Rv^ - R'Lo — 



R 



-^R'v 



2 

dE I 5 



dR 



dR 



^R'v 



'dR 



5% 
dR^ 



dR \12 

d^E 



^ 1 5£o ^ -e, dv 
vio - -Rv — --2RI.Q — 
" 2 dR dR 



+ vI.oR 



dR^ 



cos(m^ + (fio) 



again considering relation (33). 

Solving for v^2 and inserting in ( |5^ yields 



2Q 



5f2 



dtdip 



dip^ 



Vr2 



H — 1LqD.{U\ cos(mip + (fio) + U2 sin(2my? + 2(/3())) - 



2 
with 

t/i(T,-R) 
E 



= - + 3/?^ — r-2/?,^ — -^ + 3 — - — 
7? \ dR dR^ % dR % dR dR 



dE I 25 dv 16 v 5Eo\ 

+ — 2v + —R — + —R- -\ 

dR\ 3 dR 3 j:q dR) 

2 d^E dE 
- -vR — - - 2R — 

3 dR^ dT 



and 



U2(t,R) = -niR^QE^ 
A solution of (pM is 



.dE E 

m'R 



VR2 = 5l[v«2e'('^^-^^)] 

1 1 
- — Ui<psm(m<p + ^0) ^ 2 sin(2mi^ + 2i^o) 

with vr2 - vr2{t, R), and with 



(52) 



(53) 



(54) 



(55) 



(56) 



with 



^'ip2 



+ -^Wi sin(m(fi + 1^0) ~ —mU iip cos(nup + ipo) 
+ W2 ( cos(2m(^ + 2^0) +1) 

1 I dvRi _ dvR\ \ Rv d Ivr\\ 
2Q\ dT dR j 3QdR\R^) 

2 d 



3RtQD. dR 
1 / tj 



Wi(t,R) 

^ e' 

R 



dv 

6mv - (2m +\)R— + (6m - 3)7?^^ 



dR 



dR^ 



I.^^ dR ^0 oR 



1 V 

6mR^ — 
^0 



5S„ 



dE_ 

dR 



dR 
17m - 6 



+ 6m7;2 



to dR^ j 
26m - 25 dv 



-V + 



3 3 

V 

26m - 16 V ^£0^ 



■'d-R 



-R— 

3 So 5/; 

/2 8m\ 5£ 
+ - + — vTJ — - - 2(m - 1)R — 
'3 3 dR^ 5t 



and 

W2iT,R) 



1 



,dE 5E 
4R 



(58) 



(59) 



(60) 
(61) 



(62) 



Inserting ( ^6[ ) and ( p8| ) into the continuity equation ( |50| ) 
again results in a secular term of the form 



dm dcL>\ 
'^m'^dRl 



(63) 



which has to vanish. Like in the 1st order case, this can only be 
achieved by 



cj = 



(64) 



Again this leads to m^ = 1 . Hence follows that all quantities in 
2nd order, v^2, vr2 and £2, are also independent of the dynami- 
cal time f, and £2 also can be integrated directly. 

There exists a second secular term in (pO|), which is appear- 
ing already in ( |5^ and (|8|) and which is proportional to Ui<fi. 
Therefore, U 1 has to vanish. 



0) = Q.{m ± 1) 



(57) Uy=Q 



(65) 
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This leads to a differential equation for E{t, R), 



dE 
dr 



E ll dv 3 d^v 
R\2dR 2^5^ 



V C?£() 



3 R dv dl.0 
2%dR~dR 



^ dR\R^ ~6dR^ 3%^~dR)~ 3''W' ^ 

It is worth no ting that ( |66| ) agrees exactly with Eq. (58) of 
Ogilvie ( 2001 ) in the case of a pressureless viscous fluid with 
constant kinematic viscosity v. 

To summarise, according to (^4|), all quantities up to sec- 
ond order are independent of the dynamic time f. The 2nd 
order quantities depend on the azimuthal angle <f in a linear 
combination of harmonic oscillations with frequencies nvp, nvp 
and Invp, where m - +\ and in - +\. Therefore, if a non- 
axisymmetric perturbation occurs, it has to be a superposition 
of one-armed and two-armed spiral structures in second order 
This is an effect of the non-linearity of the analysis. 



4. Local stability 



The solution of (66) rules the stability or instabiUty of the vis- 
cously spreading ring. Assuming that short wavelengths be- 
come first unstable (which is supported by the numerical re- 
sults), we make the approach 



£(T,/;) = ^R[£oe'^^^-^^)] 



(67) 



where £0 is a constant. Then, (g6|) gives the dispersion relation 



with 

_\ dv 3 d'^v V d% 3 R dv d% 
Vi{t,R) - -— + - Y^'dR ^ 2%'dR~dR 



(68) 



(69) 



and 

y,(,,^).r + ^|r + !^f£ (70) 
^ ' R 6 dR 3% dR 

which are real functions. As soon as the imaginary part of cr 
becomes positive. 



^+^^>0 



(71) 



R 3 

the ring becomes unstable. This inequality implies that the 
axisymmetric ring is always unstable to perturbations of suf- 
ficiently short wavelength. Since there exists also a non- 
vanishing real part of cr, which evokes oscillatory behaviour, 
and because the instability is driven by the ki nematic viscosit y 
V, it is a type of viscous overstability (see e.g. Kley et al. 1993 ). 

Assume v - const. Then, the dispersion relation (|68|) takes 
the form 



/ V 8 V c?£o 
a- - -k\ - + ----- 

\R 3 So dR 



+ /U - - 



V dl.0 
3 Rto dR 
and the condition for instability is 
3 dto 



Rto dR 



(72) 



(73) 



Using the analytic solution (J^ of the surface density £0, the 
growth rate o",- becomes 



CT; - \k^v + 



IRlx^T* 



(t* + - Ax 



(74) 



Therefore, with the assumption of small viscous times, i.e. 
T* <K X, the growth rate is 



1,2 3 y 1 / Ro\ 

cr,- - -k v + T + — 1 

3 4/;2 6t\ r) 



(75) 



where suitable approximations of the modified Bessel func- 
tions have been applied. The condition for instability is accord- 
ingly 



4/?2 2vt\ Rl 



(76) 



With the assumption of large viscous times, i.e. t* » x, the 
growth rate becomes 



1 , 1 
/^-^67 



(77) 



Obviously, for late times the viscous ring is unconditionally 
unstable against non-axisymmetric perturbations. 

Although according to the dispersion relation the instability 
should grow faster with increasing wave number, in the numer- 
ical simulations finite wavelengths dominate the perturbation. 
There are two reasons why the instability might not occur at 
very short wavelengths. First, the relevant scales may not be 
resolved in the numerical calculations. This effect is indicated 
in Fig. |l^, which shows that the instability grows faster with 
higher spatial resolution of the numerical method. Secondly, 
the analysis will break down when the separation of dynamical 
and viscous timescale becomes invalid, i.e., below the charac- 
teristic radial scale R ~ (v/Q)'^^. 

The main limitation of the presented analysis is the omis- 
sion of pressure. In particular, the more general analysis of 
pgilvie (2001) shows that, in a disc with significant pressure, 
the eccentricity vector precesses rapidly and differentially due 
to pressure, in addition to its viscous growth or decay. When 
pressure dominates over viscosity, the afore mentioned char- 
acteristic radial scale where the analysis may break down is 
replaced by the semi-thickness of the disc. Furthermore, the ec- 
centric instability can be enhanced or suppressed by allowing 
for the kinematic viscosity to depend on the surface density, by 
introducing a bulk viscosity, by allowing for stress relaxation, 
or by including three-dimensional effects. 

5. Numerical models 

To verify the theoretical results obtained in Sects. ^ and ^ 
we performed various numerical simulations of the viscously 
evolving ring, concentrating on the case of constant kinematic 
viscosity v, = const. To be able to distinguish between numer- 
ical and physical effects, we selected two fundamentally differ- 
ent numerical methods for the simulations, and we inspected a 
wide range of parameters. 
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5.1. Using SPH 

The first method we used was smoothed particle hydrodynam- 
ics (SPH). This Lagrang i an par ticle method was introduced 
by Gingold & Monaghan ( 1977 ) and Lucy (1977) to simulate 
compressible flows with free bou ndaries. Fo r an ov erview of 
the SPH method see, for example, Monaghan ( 1992 ). 

In our simula tions we did not use the ar tificial viscosity 
approach (see e.g. Monaghan & Gingold 1983 ) usually applied 
to model viscid fluids with SPH, but we chos e a version of SPH 
that includes the en tire viscous stress tensor ( ^^lebbe et al. 1994 ; 
Riffert et al. 1995 ). This implementation allows for the proper 
treatment of viscous shear flows. Additionally, in compliance 
with the assumption of a cold disc, pressure was completely 
neglected in all the simulations. A detailed description of the 
application of t his SPH approach to the viscously evolving ring 
can be found in |Speith & Riffei^ ( |1999| ). 

The simulation presented below was performed with the 
following parameters. The initial SPH particle distribution rep- 
resents the surface density (^ of the analytic solution at the 
viscous time t* = 0.018, with Rq = I Rq, disc mass M = 
10"'" Mq, and central mass = 1 Mq. The kinematic viscos- 
ity was set to v., = 3 x 10"^/?|/s = 1.5 x lO'^'cm^/s throughout 
the whole simulation. The initial distribution consists of 40 000 
particles which are arranged symmetrically to make sure that at 
the beginning the centre of mass lay exactly at the origin. Due 
to the stochastic nature of the SPH method, no further seed per- 
turbations are required to bring potential instabilities to grow. 

Fig. |l| shows the evolution of the viscous ring during the 
SPH simulation in a sequence of plots of the surface density 
distribution. The top left plot represents the initial density dis- 
tribution at viscous time t* - 0.018. Immediately after the start 
of the simulation, symmetric structures develop, which can still 
be seen in the distribution at r* = 0.054 in the top middle panel. 
These narrow ring-like structures may be numerical artifacts 
due to the relaxation of the initial particle distribution, or they 
m ay be caused by th e radial viscous overstability discussed e.g. 
in Kley et al. (1993), but they vanish in the further course of the 
simulation. Instead, the development of a spiral pattern begins, 
in this case at the inner edge of the disc, until finally a distinct 
spiral structure is covering the whole face of the disc, as can 
clearly be seen in the distribution at r* = 0.270 in the bottom 
right panel of Fig. ^ 

To determine the nature of the spiral structure, a Fourier 
analysis of the evolving ring was performed during the simula- 
tion. Fig. H gives the time evolution of the first four azimuthal 
Fourier modes (m = 1 to m = 4) at radius R - R^). For small 
viscous times, the even modes m = 4 and especially m - 2 
are high, while the odd modes m - \ and m = 3 are low. This 
effect is due to the symmetric initial particle distribution and 
disappears eventually when the symmetry of the particle distri- 
bution is lost. Later, the (m - l)-mode increases exponentially, 
while all the other modes decrease or stay small. In particular, 
the {m - 3)-mode remains in the large noise background of 
the simulation, which is intrinsic to the SPH method, and the 
(m - 4)-modes decreases down to the noise background (as all 
inspected higher modes do). The (m - 2)-mode seems to end 
at a level slightly above the noise background, although this is 
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Fig. 3. Ratio of the amplitudes of the errors of azimuthal and 
radial velocities at r* - 0.270 for the disc shown in the bot- 
tom right panel of Fig. |l]. Plotted is the azimuthally averaged 
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'''Kepler I / 1 ^'R ^'spreading 



over radius, where and 



v's are the simulation results and v'Kepier and v'spi-eading are the ve- 
locities (^) and ( [TTl ) of the analytic solution. The curve roughly 
matches the expected value of 4 (dotted line). 
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Fig. 4. Absolute value of the phase shift between the errors of 
azimuthal and radial velocities, - VKepier and - Vspreading, 
for the disc at t* = 0.270 shown in the bottom right panel of 
Fig. |l], azimuthally averaged and plotted over radius. The curve 
roughly matches the expected value of | (dotted line). 



not clearly to determine because of the large oscillations of the 
modes which are also caused by the SPH noise level. 

The results of the Fourier analysis, i.e., the development 
of a dominant one-armed spiral structure with an additional 
weaker component of a two-armed spiral structure, are in com- 
plete agreement with the theoretical results of the first and sec- 



ond order perturbation analysis as derived in 3.4 and 3.5 



Further predictions of the perturbation analysis in Sect. ^ 
that can be tested easily, are the relations between the first or- 
der velocities vri and v^i. According to (p3|), their amplitudes 
have a ratio of 2, and they obey a phase shift of |. To ver- 
ify this, we more thoroughly studied the most evolved particle 
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Fig. 1. Evolution of the surface density of the viscous dust ring in an SPH simulation. Parameters of the disc aie Ro - I Rq, 
disc mass M - lO^'^M©, central mass - I Mq, and kinematic viscosity v, = 3 x IQ^^R^/s = 1.5 x lO'^cm^/s. The initial 
density distribution at viscous time r* = 0.018 is shown in the top left plot. The symmetric structures seen in the distribution at 
T* - 0.054 may be of numerical nature due to relaxation effects of the initial particle distribution. Later, the development of the 
spiral instability can be seen very clearly. 



distribution of the SPH simulation presented in Fig. i.e. the 
ring at viscous time t* - 0.270, whose density distribution is 
shown in the bottom right panel of Fig. |l} In Fig. || the rela- 
tion - VKepierl/|v« ' ^spreadingl IS plotted ovcr radius. Here, 
and denote the simulation results, which are averaged az- 
imuthally using the SPH formalism, and VKepier and Vspreading 
are the velocities according to the analytic solutions (g) and 
(11) of the viscous ring model. Taking - v'Kepiei ~ t^Vi 
vr - Vspreading ~ Vfii, the curve in diagram^ matches the value j 
sufficiently to satisfy relation (^). Accordingly, Fig. ^ shows 
the phase shift between - v'Kepier and vr - Vspi-eading, and again 
the plot matches the expected value of | satisfactorily. 

Finally, we want to compare the local stability analysis and 
the resulting dispersion relation established in Sect. ^ with the 
SPH simulation results. From the top left panel of Fig. ^ the 
growth rate at radius R = TJq of the first Fourier mode may 
be estimated to cr, = 1.8 x lO^^sec"' = 0.028 O(7;o), the lat- 
ter measured in units of the angular velocity (|o|) R - Rq. 
With dispersion relation (^5]), the corresponding wave number 
can be estimated to A: = 42/?^', which results in a wavelength 
of A - O.ISRq. Comparing that result with the surface density 
distributions of Fig. |l] shows good agreement. That can be seen 



in particular in Fig. ^ where the surface density distribution of 
the ring at viscous time r* - 0.270 is plotted in polar coor- 
dinates, and where additionally lines of constant kR + ip with 
k - 42/?Q ' are drawn in for the region around R x R^. The lines 
match well the slopes and the wavelengths of the spiral struc- 
tures. Note by the way that in this presentation of the surface 
density clearly can be seen that in this case the spiral structure 
consists of a leading (outer part of the disc) and a trailing (inner 
part of the disc) spiral. 

5.2. Using a finite difference metliod 
5.2.1. Numerical considerations 

In addition to the particle simulations above we performed runs 
using a finite difference method for solving the hydrodynamic 
equations. The code is based on the hydrodynamic program 
RH2D suited to study general two-dimensional systems includ- 
ing radiative transport ( ECley 1989). RH2D uses a fixed Eulerian 



grid and is a spatially second order accurate, mixed explicit 
and implicit method. Due to an operator-splitting approach, 
the method is partly centred in time and is therefore also in 
time in real terms of higher accuracy than the formal first or- 
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Fig. 2. Time evolution of the first four Fourier modes at R - Rq of the simulation shown in Fig. |l[ Because of the symmetric 
initial particle distribution, even modes are high and odd modes are low at the beginning. Later, all modes decrease except the 
(m = l)-mode which increases exponentially as expected. The (m - 3)- and (m - 4)-mode merge with the noise background 
of the simulation, while the (m = 2)-mode stays slightly above the noise level. The dashed line in the top left panel gives an 
approximated fit of the growth of the first Fourier mode to cr, = 1 .8 x lO^^sec"' - 0.028 0(/?o). 
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Fig. 5. Surface density distribution of the disc at t* = 0.270 as 
shown in the bottom right panel of Fig. |l]but plotted in polar 
coordinates (radius R over azimuth tp) instead of Cartesian co- 
ordinates. Additionally, lines of constant A:/? + with k - 42/?^' 
are drawn in around R x Rq. Note that the spiral structure con- 
sists of a leading and a trailing spiral. 



der. The advection is computed by means of the second order 
monotonic transport algorithm, introduced by van Leer (1977), 
which guarantees global conservation of mass and angular mo- 
mentum. Advection and forces are solved explicitly, while the 
viscosity is treated either explicitly or implicitly. The formal- 
ism for application to thin discs in (R - (p) geometry has been 
described in detail in [Kley| ( |1998| , |1999| ). 

To model the viscously evolving ring we work in cylindri- 
cal coordinates and solve exactly the hydrodynamic equations 
as given in Eqs. (|l]) to (^), using the full viscous stress tensor 
presented in (Q) to (^, using an extremely small pressure with 
Cs/v'Kepier ~ 10"^. The equations are solved in dimensionless 
units using an arbitrary basic unit length Rq. The unit of time to 
is chosen such that 



to = 



(78) 



i.e. in these units GM^ = 1, and the period of one Keplerian 
orbit at (the dimensionless radius) R - \ \s Pq - 2n. The only 
relevant physical constant entering the equations is the dimen- 
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sionless coefficient of the kinematic viscosity v.,, which is mea- 
sured in units of vq = ^g/fo- 

All the grid-based models presented in this paper use the 
same computational domain D represented by [^min,^max] x 

[y^min, V'max] with R^in = 0.2, R^^^ = 1.8, l^tnin = 0.0 and 

^max - 27T. The ring D is covered in the basic model by 
128 X 128 grid cells. However, to analyse resolution and numer- 
ical effects we use different griddings ranging from 64 x 128 to 
256 X 733. For our basic reference model we use for the con- 
stant dimensionless viscosity v, - 4.77 x 10"^, which is in fact 
identical to the value used in the particle simulations. 

The initial setup consists of an axisymmetric density profile 
following the analytic density Ering as given in Eq. (||). Because 
the original 5-function distribution is hard to represent numer- 
ically, we start, as above, with an initial density profile given 
here by the analytic solution at the viscous time r* - 0.016, 
i.e. S(r = finit, Jc = R/Ro) = Iring(r* = 0.016, jc). The slight dif- 
ference to the value 0.018, used in the SPH calculations, is of 
no importance for the subsequent evolution. For the standard 
viscosity 4.77 x 10"^ this refers to the initial (dimensionless) 
time finit - T*Rl/(l2vs) - 4.34fo- Stated differently, one vis- 
cous evolution time corresponds to 278 orbital periods Pq for 
the given viscosity. The total dimensionless mass M in the disc 
is normalised to unity. The initial radial velocity is set to zero 
and the angular viscosity to VKepiei - 

The code is written such that axial symmetry is preserved 
exactly for an explicit viscosity-solver, i.e. purely axisymmetric 
initial conditions remain exactly axisymmetric. Hence, the ini- 
tial density is disturbed randomly by typically 1 or 0. 1 % to sup- 
ply seed perturbations, and the subsequent longterm integration 
has to follow the evolution on a viscous timescale which is typ- 
ically about a: 100 orbital periods for the standard model. 



5.2.2. The standard model 

An analysis of the properties of the standard model is pre- 
sented. This refers to fixed physical parameter and initial con- 
ditions as outlined in the previous section. Starting from the 
disturbed axisymmetric density distribution the configuration 
evolves as presented in Fig. ^ where the density distribution 
is plotted for a model with a resolution of 256 x 256 and an 
initial density perturbation of 0.1%. During the initial phase 
axisymmetric disturbances appear which vary on short dynam- 
ical timescales. These radial oscillations may be rela t ed to the 
vi scous overstabilit y discussed for example in Kato ( 1978 ) or 
in Kley et al. ( 1993 ). The oscillations are damped however, and 
the system becomes axially symmetric again. Apart from these 
initial variations, there are no further indications of a variation 
on the dynamical timescale during the remaining evolution. At 
/ a 35 a trailing one-armed spiral density wave becomes visible 
in the inner parts of the ring (at R » QARq) which propagates 
slowly outwards. At later times a leading spiral appears in the 
outer parts of the ring which has a negative density gradient. 
This is in very good agreement with the SPH simulations pre- 
sented above, and with the analytical results which indicated 
the possibility of leading and trailing waves. The outer leading 
spiral is not as tightly wound as the trailing inner one. That is. 



w .15 



t ■ h 



Model 6 
f=52.4 
Avefaged 
at - T 




Fig. 7. Radial density distribution at = tt (crosses) and az- 
imuthally averaged (triangles). The averaged distribution fol- 
lows exactly the analytical zero order result of the diffusion 
equation Sdng, Eq. (^. Time is given in orbits at R - IRq. 

the radial wave number of the outer wave is smaller than the 
one on the inner side. This different behaviour may be caused 
by a wrapping up of the inner spiral by the differential rotation. 
At a given radius the radial separation of two wave crests (tight- 
ness of the spiral) is approximately constant with time, while 
for a given snapshot (time) there appears to be the tendency for 
the spirals to widen for larger radii. 

In Fig. ^a radial cut through the density at the fixed angle 
- n (crosses) is shown for the standard model (128 x 128 
grid points), for t - 52.4. Clearly seen are the different crests 
of the spirals. The radial wavelength A - 2n/k of the spiral 
is of the order of a few grid points only. Overlaid we plotted 
the azimuthally averaged density profile (triangles) which fol- 
lows exactly the analytical solution (solid line). This feature 
of the averaged density is a result of the modal form oc e'"^ 
of the developing spiral structure. The oscillations at the inner 
boundary are caused by the artificial inner boundary condition 
which allows for no outflow (vs = at /?min) for this particu- 
lar model. Opening the inner boundary yields smooth density 
distributions. The obtained growth rates are, however, not influ- 
enced by the exact boundary conditions, as long as the density 
there remains small with respect to the central density. 

In Fig. I the growth rates for the first 5 modes m - 1 upto 
m = 5 are displayed. The vertical axis refers to the decimal 
logarithm of the amplitude. Clearly seen is the random initial 
perturbation at the level of 1 % which is reflected in a start of 
all modes at the level of 10"^. During the initial evolution the 
rings spreads slowly, the amplitudes decline until at later times 
(f = 15-20) them = 1 mode begins to grow exponentially with 
time. From the plot and additional runs we may approximately 
determine the growth rate cr,- = lm( cr) for this m = 1 mode of 
the standard model to be 

cr, * 0.035 Q(7?= 1) . 
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Fig. 6. Grey-scale plot of the density in polar coordinates at different evolutionary times for a model with a 256 x 256 resolution. 
Times are indicated in orbital periods Pq at R - IRq. To compare with Fig. [l|: One viscous time refers to 278Po- The time 
evolution of the system proceeds similar to the particle simulations. In the initial phase axisymmetric perturbations are seen first 
(f = 17.45). They vanish later, and the system returns to axisymmetry (f = 30.54). Then, from the inner parts the instability 
develops as a trailing spiral. At the end of the simulation a leading spiral appears in the outer region of the ring. 



This is stated in units of the dimensionless time which is iden- 
tical to the Keplerian orbital frequency O at the radius R - IRq. 
To verify the numerical robustness of this result we varied sev- 
eral numerical parameter such as rotating frame, implicit vis- 
cosity, directional splitting, inner and outer boundary condi- 
tions, and found no significant variation. 

Of course, since the growth rate cr, = cr,(A;) depends on the 
azimuthal wave number k and thus on the grid size as well, we 
do expect a dependence on the numerical resolution. 

5.2.3. Varying resolution and viscosity 

To test the influence of the numerical resolution we first fixed 
the number of grid points in the radial direction (to 128) and 
then varied the number of the angular grid points (from 128 to 
370). The results in Fig. ^ indicate that the growth rate of the 
m - I mode (measured atR - IRq) is nearly independent of the 
resolution in (p and converges for large A^^ to one value. This is 
consistent with our result that to second order the growth rates 
should be independent of the azimuthal wave number m. 

In the next set of models the number of radial grid points 
was allowed to vary as well. The time evolution of the total 



radial kinetic energy is displayed in Fig. [T^ for models with 
different radial (and azimuthal) grid resolutions. For very small 
radial resolutions Nr = 64 there is very little or no growth. 
Then, for an increasing number of radial grid points the growth 
rates increase as well. This effect of faster growth for higher 
resolution is not an artifact of the numerical analysis but is 
to be expected from our analytical estimates. It was shown 
that the growth rates cr depend on the radial wave number k: 
The smaller the wavelength of the perturbation the faster the 
growth rates. But the smallest wavelength that can be resolved 
by numerical computations depends naturally on the resolution. 
Notice, that the wave crests in such a simulation are always re- 
solved by only a few grid cells (see Fig. Q)- 

To study the influence of the kinematic viscosity we var- 
ied V from the standard value to higher and smaller values. In 
Fig. |ll| results are shown for several models. A measure of the 
global deviation AL of the density from the axisymmetric struc- 
ture has been plotted, with AS defined as 



AS = max 

£> \2(r, (fi) + E(r, (p + n)) 
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Fig. 8. Time evolution of the first 5 azimuthal modes for the 
standard model measured at the radius R = 1. 
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Fig. 9. Time evolution of the logarithm of m = 1 mode (mea- 
sured at R = IRq) for models with varying resolution in the 
angular direction. The radial number of grid points is fixed to 
128. 



Fig. 10. Time evolution of the total radial kinetic energy (di- 
mensionless units) in the system for different numerical resolu- 
tion, as indicated by the labels. The physical parameter of the 
models are those of the standard model. 
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Fig. 11. Time evolution of the asymmetry AS, Eq. ([79[), for dif- 
ferent kinematic viscosities. Starting from the reference value 
(standard, v - A .11 x 10"^) we lowered and increased the value 
of V by the factors given in the key. 



where the maximum is taken over the whole computational do- 
main T). This measurement of AS is often a better tool to anal- 
yse the growth rates than just studying the m - 1 mode at some 
specific radius as done above. 

It can be seen that upon increasing the viscosity the growth 
rates are also increasing, which is again in agreement with our 
analytical estimates. The initial increase of AE for the two mod- 
els with the lowest viscosities (v, and 1 /2 v,) is related to the 
nearly ring like disturbances described before. These decline 
first, and then later the growth of the spiral disturbance sets in. 



The exact behaviour of the curves depends on the initial ran- 
dom perturbation as well. 

Finally we display the obtained growth rates for different 
viscosity coefficients in Fig. |l^ where the solid line denotes 
the results computed by the finite difference method. For larger 
viscosities the growth rates are also larger which is in agree- 
ment with Eq. (|7^. To compare exactly with the analytical re- 
sults, the wavelength of the most unstable mode needs to be 
known as well. However, the wave number is a function of ra- 
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Fig. 12. The growth rates (in units of Q at 7? = 1) versus kine- 
matic viscosity in units of the standard viscosity (4.77 x 10"^) 
for different kinematic viscosities and for both numerical meth- 
ods. 

dius, making a detailed comparison not as easy. It is found in 
general that the obtained growth rates (as given for example in 
Fig. 12) are in the right order with the first term of Eq. ( |7^ 
but are consistently lower The radial variation of the quantities 
may easily account for this slight discrepancy. 

Additionally, results of similar SPH calculations are plotted 
in Fig. |l2|, indicated by the dashed line. In all cases, the SPH re- 
sults fit the dispersion relation (^Sj) well within the accuracy of 
the SPH simulations. However, the growth rates as well as the 
associated wave numbers of the most unstable mode are lower 
than the grid-based results. Moreover, the growth rate increases 
slower with increasing viscosity. This behaviour may demon- 
strate that the relation of growth rate and wave number is de- 
pending on spatial resolution, initial condition, and on the nu- 
merical method. The discrepancy may also be caused by differ- 
ent evaluation methods. While for the SPH results the growth 
rates were measured locally from the evolution of the ampli- 
tude of the m - I Fourier mode at radius R = Rq, for the grid- 
based results they were determined globally by analysing the 
change of AS according to Eq. (^9|). This may also be another 
cause why the grid-based results fulfil the dispersion relation 
( 75 ) not exactly. 



6. Conclusion 

We have shown that the spiral instability found in various sim- 
ulations of the viscously evolving dust ring can be understood 
in terms of a secular spiral instability driven by viscosity. To 
accomplish this, we performed a perturbation analysis using 
a time-stretching transformation. It turned out that the widely 
known analytic solutions for the surface density evolution and 
the azimuthal velocity of the viscously spreading ring consist 
of the zeroth order terms of the expansion, while the solution 
for the radial velocity is given by the 1st order terms. 



As a result of the perturbation analysis we found that the 
ring may develop one-armed spiral structures to 1st order and 
one- and two-armed spiral features to 2nd order. A local stabil- 
ity analysis of an eccentricity function E provided a dispersion 
relation that led for a constant kinematic viscosity to a growth 
rate for secular spiral instabilities of 

cn^^k^v + 0(v/R^) 

in the limit kR '» I. The exact form of the growth rate depends 
on the viscosity prescription v = v{R). Our theoretical results 
could be confirmed in several simulations using two different 
numerical methods. 

As a consequence for numerical simulations of accretion 
discs, the spiral instability has to be taken into account when 
using the spreading ring for instance as test problem for code 
development. Otherwise, occasionally emerging spiral features 
may be mistaken for numerical instabilities of the algorithms 
used for the calculation of the viscous forces. 

The physical implications of our results are not as obvi- 
ous. Because the kinematic viscosity v is held axisymmetric 
throughout the whole stability analysis, the discovered non- 
axisymmetric instabilities may be of importance mainly in ac- 
cretion discs where the viscosity is determined uniformly by 
external effects like irradiation. 
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